{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unit convenience functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, REBOUND offers simple functionality for converting units.  One implicitly sets the units for the simulation through the values used for the initial conditions, but one has to set the appropriate value for the gravitational constant `G`, and sometimes it is convenient to get the output in different units.\n",
    "\n",
    "The default value for `G` is 1, so one can:\n",
    "\n",
    "a) use units for the initial conditions where `G=1` (e.g., AU, $M_\\odot$, yr/$2\\pi$)\n",
    "\n",
    "b) set `G` manually to the value appropriate for the adopted initial conditions, e.g., to use SI units,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rebound\n",
    "import math\n",
    "sim = rebound.Simulation()\n",
    "sim.G = 6.674e-11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c) set rebound.units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "G = 39.476926421373.\n"
     ]
    }
   ],
   "source": [
    "sim.units = ('yr', 'AU', 'Msun')\n",
    "print(\"G = {0}.\".format(sim.G))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you set the units, REBOUND converts `G` to the appropriate value for the units passed (must pass exactly 3 units for mass length and time, but they can be in any order).  Note that if you are interested in high precision, you have to be quite particular about the exact units.  \n",
    "\n",
    "As an aside, the reason why `G` differs from $4\\pi^2 \\approx 39.47841760435743$ is mostly that we follow the convention of defining a \"year\" as 365.25 days (a Julian year), whereas the Earth's sidereal orbital period is closer to 365.256 days (and at even finer level, Venus and Mercury modify the orbital period).  `G` would only equal $4\\pi^2$ in units where a \"year\" was exactly equal to one orbital period at $1 AU$ around a $1 M_\\odot$ star.\n",
    "\n",
    "**Adding particles**\n",
    "\n",
    "If you use `sim.units` at all, you need to set the units before adding any particles.  You can then add particles in any of the ways described in [WHFast.ipynb](../WHFast).  You can also add particles drawing from the horizons database (see [Churyumov-Gerasimenko.ipynb](../Churyumov-Gerasimenko)).  If you don't set the units ahead of time, HORIZONS will return initial conditions in units of AU, $M_\\odot$ and yrs/$2\\pi$, such that `G=1`.  \n",
    "\n",
    "Above we switched to units of AU, $M_\\odot$ and yrs, so when we add Earth:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Searching NASA Horizons for 'Earth'... Found: Target body name: Earth-Moon Barycenter (3).\n",
      "v = 6.370350510017522\n"
     ]
    }
   ],
   "source": [
    "sim.add('Earth')\n",
    "ps = sim.particles\n",
    "import math\n",
    "print(\"v = {0}\".format(math.sqrt(ps[0].vx**2 + ps[0].vy**2 + ps[0].vz**2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we see that the velocity is correctly set to approximately $2\\pi$ AU/yr.\n",
    "\n",
    "If you'd like to enter the initial conditions in one set of units, and then use a different set for the simulation, you can use the sim.convert_particle_units function, which converts both the initial conditions and `G`.  Since we added Earth above, we restart with a new `Simulation` instance; otherwise we'll get an error saying that we can't set the units with particles already loaded:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------\n",
      "REBOUND version:     \t3.15.0\n",
      "REBOUND built on:    \tFeb  8 2021 15:12:45\n",
      "Number of particles: \t2\n",
      "Selected integrator: \tias15\n",
      "Simulation time:     \t0.0000000000000000e+00\n",
      "Current timestep:    \t0.001000\n",
      "---------------------------------\n",
      "<rebound.particle.Particle object at 0x7f9ea41b57c0, m=1.0007667100237814 x=0.0 y=0.0 z=0.0 vx=0.0 vy=0.0 vz=0.0>\n",
      "<rebound.particle.Particle object at 0x7f9ea41b55c0, m=3.002300130071344e-06 x=1.0026880683402668 y=0.0 z=0.0 vx=0.0 vy=6.277053341010209 vz=6.277053341010209>\n",
      "---------------------------------\n"
     ]
    }
   ],
   "source": [
    "sim = rebound.Simulation()\n",
    "sim.units = ('m', 's', 'kg')\n",
    "sim.add(m=1.99e30)\n",
    "sim.add(m=5.97e24,a=1.5e11)\n",
    "\n",
    "sim.convert_particle_units('AU', 'yr', 'Msun')\n",
    "sim.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first set the units to SI, added (approximate values for) the Sun and Earth in these units, and switched to AU, yr, $M_\\odot$.  You can see that the particle states were converted correctly--the Sun has a mass of about 1, and the Earth has a distance of about 1.\n",
    "\n",
    "Note that when you pass orbital elements to sim.add, you *must* make sure `G` is set correctly ahead of time (through either 3 of the methods above), since it will use the value of `sim.G` to generate the velocities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "G = 1.0\n",
      "---------------------------------\n",
      "REBOUND version:     \t3.15.0\n",
      "REBOUND built on:    \tFeb  8 2021 15:12:45\n",
      "Number of particles: \t2\n",
      "Selected integrator: \tias15\n",
      "Simulation time:     \t0.0000000000000000e+00\n",
      "Current timestep:    \t0.001000\n",
      "---------------------------------\n",
      "<rebound.particle.Particle object at 0x7f9ea41b5540, m=1.99e+30 x=0.0 y=0.0 z=0.0 vx=0.0 vy=0.0 vz=0.0>\n",
      "<rebound.particle.Particle object at 0x7f9ea41b5d40, m=5.97e+24 x=150000000000.0 y=0.0 z=0.0 vx=0.0 vy=3642349031.417317 vz=3642349031.417317>\n",
      "---------------------------------\n"
     ]
    }
   ],
   "source": [
    "sim = rebound.Simulation()\n",
    "print(\"G = {0}\".format(sim.G))\n",
    "sim.add(m=1.99e30)\n",
    "sim.add(m=5.97e24,a=1.5e11)\n",
    "sim.status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The orbital speed of Earth is $\\sim 3\\times 10^4$ m/s, but since we didn't correctly set `G` ahead of time, we get $\\sim 3\\times 10^9$ m/s, so the Earth would fly off the Sun in this simulation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
